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We present a theoretical framework to analyze the dynamics of gene expression with stochastic 
bursts. Beginning with an individual-based model which fully accounts for the messenger RNA 
(mRNA) and protein populations, we propose a novel expansion of the master equation for the joint 
process. The resulting coarse-grained model reduces the dimensionality of the system, describing 
only the protein population while fully accounting for the effects of discrete and fluctuating mRNA 
population. Closed form expressions for the stationary distribution of the protein population and 
mean first-passage times of the coarse-grained model are derived and large-scale Monte Carlo sim¬ 
ulations show that the analysis accurately describes the individual-based process accounting for 
mRNA population, in contrast to the failure of commonly proposed diffusion-type models. 
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Intrinsic noise originating from the discreteness of in¬ 
teracting particles plays an important role in genetic ex¬ 
pression: it diversifies the distribution of protein popula¬ 
tion, promotes transition between different cellular phe¬ 
notypes on a population level, and in turn enhances or¬ 
ganisms’ ability to adapt to changing environments with¬ 
out the need of genetic mutation [1]. There are two pri¬ 
mary sources of intrinsic noise in the context of gene ex¬ 
pression: transcriptional noise from the stochastic tran¬ 
sition between active and repressed states of DNA tran¬ 
scription, and translational noise from the relatively fast 
action of mRNA to produce the proteins [1, 2]. Both 
steps result in bursts of protein production which are ex¬ 
perimentally observed [3, 4]. 

Many stochastic models have been proposed to model 
gene circuits [5-14] but only a few studies quantitatively 
account for the effects of bursting noise [11, 12, 15, 16]. 
To our knowledge, current theoretical investigation of the 
dynamical properties of such bursting processes is limited 
to stationary properties of the protein distribution on the 
population level [15, 16]. 

This Letter presents a new mathematical framework to 
analyze bursting noise in gene expression. Starting from 
an individual-based model including both mRNA and 
protein populations we construct a novel coarse-grained 
process describing only the protein population dynamics 
that fully accounts for the discreteness effects and fluc¬ 
tuations in the mRNA population. When the mRNA 
degrades at a much shorter time scale, the approxi¬ 
mating process converges to currently proposed bursting 
models [15, 16]. In our process-based framework, mean 
first-passage times in a autoregulated gene circuit with 
stochastic bursts can be formulated and calculated. 
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FIG. 1. Schematic diagrams of (a) the individual-based model 
and (b) the piecewise deterministic Markov process model (3). 
(c) Hill functions with Hill coeffcients n = 1,2, 4. 


We present analytic solutions along with computa¬ 
tional verification from large-scale Monte-Carlo simula¬ 
tions. A key conclusion is that the conventional diffusion 
approximation of the master equation fails to accurately 
estimate switching times of the individual-based model. 

A simple individual-based model of autoregulated gene 
expression including both the mRNA and protein pop¬ 
ulations contains four reaction steps [2, 15] as summa¬ 
rized in Fig. 1(a): synthesis of mRNA’s (transcription, 

(f> — P - \ mRNA), production of the protein (trans¬ 
lation, mRNA 7B °> mRNA + Protein), and degrada¬ 
tion of both the mRNA’s and proteins (mRNA </> 
and Protein -^>). In the first step IVp is the ran¬ 
dom number of available proteins and in this autoreg¬ 
ulated genetic circuit, the population of proteins regu¬ 
lates the transcription rate. The Hill function Hq(X) := 
ro + riX n /(K n + X n ) with the Hill coefficient n approxi¬ 
mates the autoregulated transcription rate when the gene 
switches between on and off states on a much shorter time 
scale [2]. 
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We refer to the process in Fig. 1(a) as the individual- 
based (IB) model. Although the IB model provides a 
detail description of both the mRNA and protein popu¬ 
lations, it is generally difficult to analyze theoretically 
except for linear cases [17, 18]. Single-species mod¬ 
els describing only the protein populations are often 
adopted, especially for more complicated genetic circuits 
[9, 10, 13, 14]. However, fluctuations in the mRNA pop¬ 
ulation are an important dynamical factor [ 11 , 12 ] and 
our objective is to construct a coarse-grained model de¬ 
scribing only the protein population accounting for con¬ 
tributions from fluctuation in the mRNA population. 

Generally mRNA’s degrade much faster than proteins. 
In the model organism Escherichia coli for example, the 
mean lifetime of the mRNA is about 2 min while protein 
lifetimes are 45 ~ 60 min [15]. As a consequence, a large 
number of proteins is produced in a relatively short pe¬ 
riod of time—a phenomenon termed translational burst¬ 
ing. In addition, due to small system size (the volume 
of E. coli are ~ 10 - 18 m 3 ), the onset of the transcription 
and the lifetime of the synthesized mRNA are observed 
in a stochastic manner [19]. 

Motivated by the observation of translational bursting, 
we propose a novel expansion to approximate the master 
equation of the IB process in Fig. 1(a). First, we no¬ 
tice that in the IB model, for any given mRNA number 
m, the protein population Np(t) is a birth-death pro¬ 
cesses with constant birth rate and constant per 

capita death rate 70 . Therefore, it is convenient to ex¬ 
pand the process describing the protein dynamics con¬ 
ditioning on the mRNA population: each “state” of the 
system is labeled by the mRNA number m. The transi¬ 
tion rate from state m to m + 1 mRNA molecules is the 
autoregulated transcription rate Hq(Np), and the tran¬ 
sition rate from state m + 1 to m mRNA molecules is 
the mRNA degradation rate 7 . Within each state of the 
system we perform a Kramers-Moyal expansion of the 
birth-death process [ 20 , 21 ] with respect to the system 
size K 1. In the lowest order approximation only the 
advection terms describing the mean-field dynamics are 
retained [22], Formally letting the protein concentration 
be x := Np/K > 0 and the number of mRNA molecules 
be m € { 0 , 1 , 2 ,...}, in each state the protein density 
evolves according to the deterministic equation 

x(t) = mj-j^ - 7 0 x, ( 1 ) 

with transition rates between different states 

H(x) '■y . . 

m - > m + 1 and m —> m — 1 ( 2 ) 

where H{x) := Hg(Kx) is the scaled Hill-function. 

Next we note that the mean lifetime of the mRNA is 
0 ( 1 / 7 ) and in the fast-degrading mRNA limit 7 1, 

most of the time the system has either m = 0 orm=l. 
We therefore neglect states m > 2 and formulate a 



FIG. 2. Sample paths of the models. Dotted green line de¬ 
notes 165 molecules which separates low and high protein 
abundance mode. 


closed forward equation for p m (x,t), the joint probabil¬ 
ity density that the system presents m G {0,1} mRNA 
molecules and protein density x at time t: 


d_ f pi{x,t) \ t ( PiOM) \ 
dt\Po(x,t)J \Po(x,t)J' 

where the forward operator [23] is defined to be 

T t , = ( -7 - d x (76 - 70 Z) H(x) 

V 7 -H(x) + jod x x 


( 3 ) 


( 4 ) 


where we have defined b := Bq/K to be a dimensionless 
parameter characterizing the strength of the bursts. We 
shall refer to (3) as the piecewise deterministic Markov 
process (PDMP; schematic diagram Fig. 1(b)) and re¬ 
mark that the process in x alone is non-Markovian. 

To proceed with our analytic investigation, an in¬ 
finitely fast-degrading mRNA limit 7 —>- 00 is taken. 
Although in such a limit the mean duration when the 
system stays in m = 1 state is I /7 —>• 0 , the protein con¬ 
centration in to = 1 state increases with a rate 67 —» 00 , 
preserving exponentially distributed random burst with 
an average burst strength b. In this limit the process 
stays in m = 0 state almost surely (i.e. p±(x,t) —0 Vt), 
and the probability distribution po satisfies a closed and 
second-order differential equation 


(1 + bd x ) d t po = ~dx [-x + bH{x) - bj 0 d x x\p 0 . (5) 


The stationary probability distribution is obtained by di¬ 
rect integration: 


. , M -x x H(£)dn 

Psta.t{x) = - exp { — + / -— } ( 6 ) 

7o£ l 0 J 7o4 J 

where M is the normalization factor. Substituting the 
explicit form of the Hill function we find the analytic 
expression for the stationary distribution 

Pst a t{x) = — e~^x^~ 1 (x n + 1 ”)^ . (7) 

7o 

We remark that in the limit 7 —> 00 the PDMP model 
reduces to the bursting model described by a continuous 
master equation, and that this result confirms [16]. 
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FIG. 3. Top panel: Drift of the mean-field dynamics x = 
bH ( x )— 70 x showing a single fixed point. Bottom panel: Sta¬ 
tionary probability distributions of individual-based model, 
piecewise deterministic Markov process (PDMP), and the dif¬ 
fusion approximations (DA). Solid lines are analytic solutions. 
Discrete markers represent numerically measured probability 
distributions from Monte Carlo simulations. 


In some parameter regimes the stationary distribution 
(7) exhibits bi-stability [24] and can be adopted to model 
a biological switch [11, 16]. Our formulation (3) can be 
used to derive the mean switching time (MST) between 
two modes of gene expression in a straightforward way 
[25]. We begin by deriving the mean first-exit time to 
leave a domain (xi,X 2 ) where 0 < x\ < X 2 < oo. 

If the initial protein concentration is x £ (xi,x 2 ) and 
the initial number of mRNA is m , then the mean time 
to exit the domain T m (x) satisfies the inhomogeneous 
equation [ 20 ] 



= L 


Ti{x,t) \ 
Tq(x, t) ) ’ 


( 8 ) 


where the generator L is the adjoint of the forward op- 

I 


erator in (4), 


L = 


-7 + (7 b - 7 qx) d x 


7 


-H(x) - 'yoxdx 


H(x) 

with boundary conditions 

Ti(x 2 ) = 0, T 0 (x 1 ) = 0. 


(9) 


( 10 ) 


The physical meaning of the boundary conditions is clear: 
when the system starts with the state m = 1 —a state 
with fast production of proteins—at upper boundary x 2 , 
and when the state m = 0 —a state with only degrading 
proteins—at lower boundary xi, immediately the flow 
leaves the domain (xi,x 2 ). 

Taking the limit 7 —>• 00 we deduce a closed second- 
order differential equation for To, 


-T" 


H 1 H / x y 

7 “ b + 7 Iff/ 


T' — 
1 0 ~ 


1 

&7 0 x 


iff 

70 xff 


( 11 ) 


where prime denote derivative with respect to x. The 
boundary conditions for ( 11 ) follow from ( 10 ): 


T 0 (xi) = 0,1 = H (x 2 ) T 0 (x 2 ) + J 0 x 2 Tq (x 2 ) • (12) 

We remark that while formally deriving the backward 
equations of the bursting models [15, 16] considering only 
the protein concentration is possible, imposing the cor¬ 
rect boundary conditions ( 12 ) is not trivial. 

The solution (derived in the Supplemental Material) is 

T 0 (x) = C f e~ M{y) dy+ [ e~ M ^V(y)dy, (13) 

J X l J X l 


where the auxiliary functions M (x) and V (x) and the 
constant C are 


M(x) 

V(x) 



' H(y) 

. f(y ) 

\hoy 


_ £ d y 

b dy\ n H(y) 
+ -1 dH(y) \ 

7 o yH(y) dy J 



(14) 

(15) 


C = 


rx 2 

-V(x 2 )e~ M ^f(x 2 ) - H[x 2 ) / V{y)e~ M ^dy + 1 

J X-i 


f{x 2 )e 


-M(x 2 ) 


H(x 2 ) 


-M(y) , 


(16) 


This solution is a generalization of results in [26, 27]. 

When the system exhibits bi-modality, the mean 
switching times between two modes of protein expres¬ 
sion can be obtained by taking appropriate limits of (13). 
First, we define a critical density x c separating the low 
and high protein abundance modes, then take xi —> 0 and 
x 2 —> x c for the low mode, and X\ —> x c and x 2 —> 00 
for the high mode. Careful analysis is needed because 


(11) is singular at x = 0 (and is presented in the Supple¬ 
mental Material). The analytic expressions for the mean 
switching times (MSTs) are 


T\ow— ^high — 


-Ihigh—flow — 


e- M ^V(y)dy + C 2 , 

e ~M(y) J y-V(oo)]dy, 


(17) 

(18) 
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where the constant C 2 is 

tf(Zc) Jo 

We now turn to the diffusion approximation (DA) of 
the IB process. To our knowledge there is no standard 
way to derive DA models for general bursting kernels. 
In the Supplemental Material we present the straightfor¬ 
ward Kramers-Moyal expansion [20, 21] of the master 
equation of the IB process in the limit 7 —> 00 yielding 
the Ito stochastic differential equation 

dX t = [bH(X t ) - 7 o W] dt + V Fb°-H (*) dW t (20) 

where X t is the random population density of the pro¬ 
teins, W t is the standard Wiener process and the scal¬ 
ing factor T = 2. An alternative and phenomenolog¬ 
ical construction the diffusion approximation is to in¬ 
sert the mean and variance of the bursting kernel in 
the individual-based process (see Supplemental Material) 
which yields (20) with the scaling factor T = 1. To avoid 
leaking probabilities to negative densities, we put a reflec¬ 
tive boundary at the origin x = 0. Analytic expressions 
for the stationary distribution and the mean switching 
times of the diffusion equation are derived by standard 
analysis [20] and presented in the Supplemental Material. 

We performed numerical simulations to measure the 
stationary distributions and the mean switching times 
(MST) in all three models to verify the theoretical anal¬ 
ysis. For the IB model, exact sample paths are generated 
by standard continuous time Markov chain simulations 
[28]. For the PDMP model, kinetic Monte Carlo simu¬ 
lations can be constructed by generating exact random 
waiting times to the next transition events [29]. In the 
limit 7 —>• 00 , we adopted a previously proposed algo¬ 
rithm [30]. For the diffusion approximations we construct 
a standard Euler-Maruyama integrator of (20). 

The parameters were chosen to be in a biologically rel¬ 
evant regime [15, 16, 31]: K = 200, n = 4, B Q = 40, 
r 0 = 2, n = 10, 7 = 30, while 70 := 1 is chosen to 
normalize the unit of the time by a natural cell cycle. 

Fig. 3 presents the stationary probability distributions 
of the IB, PDMP, and DA models. Note that the low¬ 
mode is noise induced and does not exist in the mean-field 
dynamics (top panel of Fig 3). While the PDMP model 
captures the stationary distribution of the IB model ex¬ 
tremely well, directly expanding the IB stochastic burst¬ 
ing model by Kramers-Moyal expansion (DA with T = 2) 
qualitatively captures the stationary distribution, and 
the phenomenological DA model with T = 1 failed to 
capture the stability of the low mode. 

Fig. 4 presents the MST between low and high protein- 
abundance modes in all three models. Again, the PDMP 
model well estimates the mean switching times of the IB 
model, and both the DA models fail by a large amount. 



FIG. 4. Mean switching times (MST) of the individual-based 
model, piecewise deterministic Markov process (PDMP), and 
the diffusion approximations (DA). Solid lines are analytic 
predictions, and discrete markers are measured from Monte 
Carlo simulations. Left panel: Ti ow _>high; Right panel: 
Thigh-now x c ~ 1Q5/K. 

When the state is initially below x c , both the DA mod¬ 
els under-estimate the transition time because the burst¬ 
ing kernels of the DA model have a thinner (Gaussian) 
tail compared to to the geometric (for the derivation see 
Supplemental Material) bursting kernel of the IB model. 
When the initial state is above x c , the DA model with 
r = 1 over-estimate the MST because the the approxima¬ 
tion does not capture the high probabilities of low-density 
bursts, and the DA model with T = 2 underestimate the 
MST because the approximation fail to capture that the 
bursting kernel is always positive. 

The PDMP approximation works well even for mod¬ 
els with a strong noise strength. In our example, the 
low-mode is of an order of 100 protein molecules, and 
the noise strength (per each burst) is of order 40 pro¬ 
tein molecules. In addition, the PDMP approximation 
performs well even though an infinitely-fast degrading 
mRNA limit 7 —> 00 is taken and consequently almost 
surely there is no mRNA presented in the system. Mean¬ 
while we observe an average 0.3188 mRNA in the station¬ 
ary state of the IB model. 

The PDMP model can be easily generalized. For ex¬ 
ample, finite population and lifetime of mRNA can be 
considered by generalize (3) to include p m with m £ 
{0,1,2,...}. The transcriptional bursting can be in¬ 
cluded by considering multiple stages of the gene. Higher 
dimensional genetic circuit can be investigated by includ¬ 
ing more states of the system [32]. These generalizations 
merit future investigations. 

We conclude that bursting originating from the dis¬ 
creteness of the fast-living mRNA molecules and the 
stochastic transcription events is the dominating noise 
in individual-based autoregulated gene expression model. 
Diffusion approximations are no longer adequate to ana¬ 
lyze the dynamical properties of bursting systems while 
the novel expansion described here faithfully captures the 
dynamical properties of the individual-based model in a 
biologically realistic parameter regime and serves as a 
new analytic tool to investigate more complex models 
with bursting noise. 








5 


Acknowledgement. YTL was supported by the visi¬ 
tor program at MPI-PKS and EPSRC (grant reference 
EP/K037145/1). CRD was supported in part by US- 
NSF Awards PHY-1205219 and DMS-1515161, and as 
Simons (Foundation) Fellow in Theoretical Physics. 


* yenting@umich.edu 

* doering@umich.edu 

[1] M. Kaern, T. C. Elston, W. J. Blake, and J. J. Collins, 
Nature Rev. Genet. 6, 451 (2005). 

[2] A. M. Walczak, A. Mugler, and C. H. Wiggins, Methods 
in molecular biology 880 , 273 (2011). 

[3] E. M. Ozbudak, M. Thattai, I. Kurtser, A. D. Grossman, 
and A. van Oudenaarden, Nature Genet. 31 , 69 (2002). 

[4] W. J. Blake, M. Kaern, C. R. Cantor, and J. J. Collins, 
Nature 422 , 633 (2003). 

[5] T. B. Kepler and T. C. Elston, Biophys. J 81, 3116 

( 2001 ). 

[6] J. E. M. Hornos, D. Schultz, G. C. P. Innocentini, 
J. Wang, A. M. Walczak, J. N. Onuchic, and P. G. 
Wolynes, Phys. Rev. E 72 , 051907 (2005). 

[7] A. M. Walczak, M. Sasai, and P. G. Wolynes, Biophys. 
J. 88, 828 (2005). 

[8] P. B. Warren and P. R. ten Wolde, J. Phys. Chem. B 
109 , 6812 (2005). 

[9] J. Wang, L. Xu., E. Wang, and S. Huang, Biophys. J. 
99 , 29 (2010). 

[10] J. Wang, K. Zhang, L. Xu, and E. Wang, Proc. Natl. 
Acad. Sci. U.S.A. 108 , 8257 (2010). 

[11] M. Assaf, E. Roberts, and Z. Luthey-Schulten, Physical 
Review Letters 106 , 248102 (2011), pRL. 

[12] M. Strasser, F. J. Theis, and C. Marr, Biophys. J 102 , 
19 (2012). 

[13] J. X. Zhou, M. D. S. Aliyu, E. Aurell, and S. Huang, J. 
R. Soc. Interface 9, 3539 (2012). 


[14] M. Lu, J. Onuchic, and E. Ben-Jacob, Phys. Rev. Lett. 
113 , 078102 (2014). 

[15] M. Thattai and A. van Oudenaarden, Proc. Natl. Acad. 
Sci. U.S.A. 98, 8614 (2001). 

[16] N. Friedman, L. Cai, and X. S. Xie, Phys. Rev. Lett. 97, 
168302 (2006). 

[17] V. Shahrezaei and P. S. Swain, Proc. Natl. Acad. Sci. 
U.S.A. 105 , 17256 (2008). 

[18] N. Kumar, T. Platini, and R. V. Kulkarni, Phys. Rev. 
Lett 113 , 268105 (2014). 

[19] L. Cai, N. Friedman, and X. S. Xie, Nature 440 , 358 
(2006). 

[20] N. G. van Kampen, “Stochastic processes in physics and 
chemistry,” (North-Holland, Amsterdam, 2007). 

[21] T. G. Kurtz, J. Appl. Probab. 8, 344 (1971). 

[22] T. G. Kurtz, J. Appl. Probab. 7 , 49 (1970). 

[23] Notation: the differential operator d x acts as a total dif¬ 
ferential including the probability distributions p m (x,t) 
outside the matrix. 

[24] W. Horsthemke and R. Lefever, “Noise induced transi¬ 
tions: Theory and applications in physics, chemistry, and 
biology,” (Springer, Berlin-Heidelberg, 2006). 

[25] C. R. Doering, Phys. Rev. A 34 , 2564 (1986). 

[26] J. Masoliver, K. Lindenberg, and B. J. West, Phys. Rev. 
A 34, 2351 (1986). 

[27] C. R. Doering, Phys. Rev. A 35 , 3166 (1987). 

[28] R. Schwartz, Biological modeling and simulation : a sur¬ 
vey of practical models, algorithms, and numerical meth¬ 
ods (MIT Press, Cambridge, Mass., 2008). 

[29] The waiting times can be generated by deriving the sur¬ 
vival functions of each state [26] and utilizing the inverse 
transform sampling. 

[30] P. Bokes, J. King, A. T. A. Wood, and M. Loose, Bull. 
Math. Biol. 75 , 351 (2013). 

[31] Y. Taniguchi, P. J. Choi, G. W. Li, H. Y. Chen, M. Babu, 
J. Hearn, A. Emili, and X. S. Xie, Science 329, 533 
( 2011 ). 

[32] Y. T. Lin and T. Galla, arXiv:1508.00608 [q-bio.MN] 
(2015), submitted. 



Supplemental Material of 
Gene expression dynamics with stochastic bursts: 
exact results for a coarse-grained model 

Yen Ting Lin* 

Theoretical Physics Division , School of Physics and Astronomy, 

The University of Manchester , UK and 
Max Planck Institute for the Physics of Complex Systems, Germany 

Charles R. Doeringl' 

Department of Mathematics, University of Michigan, USA 
Department of Physics, University of Michigan, USA and 
Center for the Study of Complex Systems, University of Michigan, USA 


* yenting@umich.edu 
t doering@umich.edu 


1 



I. INDIVIDUAL-BASED MODEL INCLUDING THE MRNA POPULATION 


In this section, we define the individual-based model including the mRNA populations. The 
kinetic scheme of the autoregulated process is [1, 2] 


(j)—> 1 x mRNA 
mRNA —> mRNA + 1 x Protein 
mRNA —> (j) 

Protein —> <fi 


with a rate H 0 (N Plotein ), 
with a rate 7B0) 
with a rate 7, 
with a rate 70, 


(1) 


where the Hill function is defined to be 


H q {x) := r 0 + n 


k n + x" 


( 2 ) 


We will denote the (random) population of mRNA and protein by N m and N P respectively. 
Define the joint probability distribution of the system at time t to be 


Pm,n{t) '■= P {A^ m = m, N P = 7i at time t} . 
The master equation of process ( 1 ) is 


( 3 ) 


P m ,n = - [H 0 (n ) + 7R0TO + 7m + 7 0 n] P m ,„ ( 4 ) 

+ 7 (m + 1) P m+ i,n + H 0 (n)P m - i, n + 7&mP mjra _i +70 (n+ 1) P m ,„+i. ( 5 ) 

Continuous time Markov chain simulations [ 3 ] are constructed to generate exact sample paths. 


II. DERIVING THE PDMP IN THE LARGE POPULATION LIMIT 


In the fast degrading mRNA limit ( 7/70 1), the system only presents only 1 or 0 mRNA 

in a majority portion of the time. Our proposed approximation is to consider the process (1) 
conditioning on whether or not the system presents an mRNA, and truncate the probabilities 
associated with mRNA number greater than 1: 

P m ,n = 0, Vm > 1, ( 6 a) 

Pl.n = - [jB q + 7 + 7o7i] Pi,n + H Q (n)P 0 ,n + 7P 0 Pi,„_i + 7o {ji + 1) Pi,„+i, ( 6 b) 

Po.n = - [H 0 (n) + 7 o 7 l] P 0 ,„ + 7Pl,n + 7o (n + 1 ) Pl,n+ 1 - ( 6 c) 

Next, for each of the master equations of the protein number n conditioning on the mRNA 
number m, we perform the conventional Kramers-Moyal expansion[4]. Denote a typical popu¬ 
lation scale of the protein by Nq 1. Note that in the autocorrelated circuit, it is convenient 
to choose Nq = K. In the continuum limit, the population density is defined by x := ti/A', 
and the mean “burst” size is defined is b := Bq/Nq. The evolution of the probability distri¬ 
butions po(x,t) := Po, n (t)/K and p\(x,t ) := P\^ n {t)/K is well-approximated by two coupled 
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Fokker-Planck equations [5] 


d t p+ = - 7 p+ + H(x)p- + 


d x hox - 76 ) + — d 2 x ( 70 a; + 76 ) 


P+, 


5iP- = - H(x)p- + 7p+ + 70 ( + — dlx ) p-. 


The coupled Fokker-Planck equation can be expressed in a compact matrix form: 


(7a) 

(7b) 


d t 




( 8 ) 


with 


Lt := 


-7 + dx (7ox - 7 b) + jk d 2 x (70X + 7 b) 
7 


H \ 

-H + l0 d x x + 2Nlcdlx ) 


(9) 


We again remind the reader that the differential operators d x and df. act on p 0 ,i too. 

It should be clear that the the discrete population of the proteins causes the demographic 
stochasticity, which is described by those terms with a prefactor 1/K. We further propose to 
take the limit K —» 00 [5] and leave only the advection terms in (7) to consider exclusively the 
bursting noise, which is a result of randomly production and degradation events of the mRNA. 
In such a limit, the process becomes a piecewise deterministic Markov process: in each state of 
m = 0 or m = 1, the process is deterministic but the switching between the states is Markovian. 
We emphasize that, in such a limit, the demographic noise which comes from the discreteness of 
the protein population does not exist—condition on a m state, the concentration of the protein 
on its own is always evolving in a deterministic fashion. 

We notice that the duration of m = 1 state is of order O (I/ 7 ), but the resident time of the 
to = 0 state does not depend on 7 - it is of of order O ( 7 0 ). As a consequence, pi scales O (I/ 7 ) 
and as 7 —> 00 . In fact, it can be rigorously proved that p\ —> 0 as 7 —» 00 , and jp± in (7b) can 
be eliminated by (7a): 


1 + bd x 




[H(x)p-]+ 7 o (1 + bd x ) 


—d 2 z 
2 K x 



III. INDIVIDUAL-BASED BURSTING MODEL 


Back to process ( 1 ), when 7 70 and -ff(-/Vp ro tein), there is a time-scale separation and the 

mRNA’s degrade at a very rapid rate. As a consequence, when one mRNA is formed, almost 
surely the next happening events before its final degradation are the even more rapid productions 
of proteins. 

Due to the time-scale separation, the production of other mRNA’s is and the protein degrada¬ 
tion are negligible in one mRNA’s lifetime. In such a limit, the distribution of the total number of 
proteins an mRNA could ever synthesized before its final degradation can be computed. Define 
the total number of proteins an mRNA could ever synthesize to be Vs, a non-negative random 
variable. Because the mRNA has only have two choices—either to degrade or to produce a 
protein—at any time before the final degradation, the probability that the mRNA produce a 
protein is -Bo/(l + Bq) from reading off the ratio of the rates in the process (1). Therefore, the 
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distribution of Ns is a geometric distribution 


>{7V s =n} = 


Bn 


l + B 0 


1 


1 + Br 


Ns + 1 ~ Geom 


Bn 


1 + Bn 


( 11 ) 


As a consequence, process (1) in the limit 7 —> 00 can be re-formulated to neglect the mRNA 
population 

<t>—> Ns x Protein with a rate Hq (-ZVp), 

Protein— cj) with a rate 70 , (12) 

Ns + 1 ~ Geom (B 0 / (1 + B 0 )). 


We remind the reader that the parameter Bq is the mean number of the proteins an mRNA can 
produce. We shall refer to model (12) as the “individual-based bursting model”. 

Let P n to be the probability when the system has exactly n proteins. The master equation of 
process ( 12 ) can be derived 


A = - 


n / r> \ r 

[Ho(n) + 70 n] Pn+Yl H o{m) I ° ) 




IV. EQUIVALENCE BETWEEN THE PDMP AND CONTINUOUS STATE BURST¬ 
ING MODELS 


We now apply Kramers-Moyal system-size expansion is performed only to the degradation 
dynamics in (13). The expansion of (13) yields 


d t p(x, t) 


lod x x + d 2 x x 


p{x,t) + 



y)H{y)p{y)dy, 


(14) 


where p(x,t) := P n (t)/Nn is the continuum probability distribution, x := n/Nn is the popu¬ 
lation density of the protein, and W(x — y) is a kernel of the bursting process, defined by the 
approximating the discrete by the trapezoid rule: 



y)f{y)dy 


-/m + 5 


1 ( m , m 


+ 


/ 0 1 /K + ^ 


1 + bK l + bK / 

-WV n (*-»)log(l+^ ff ) f( y ) dy ' 


In the infinity population limit K —> 00 , (14) reduces to 


(15) 


f x e— 

d t p{x, t) = d x [xp(x,t)] + J —-— H(y)p(y)dy - H(x)p(x), 


(16) 


which is exactly the continuous master equation in Friedman et al. [ 6 ] 

It is straightforward to establish the equivalence of the piecewise deterministic Markov process 
[i.e., (7) in the limit K —> 00 ] and the continuous-state bursting model (12): acting an operator 
1 + bd x to (16) yields (10) as K — > 00 . 
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V. DERIVING THE DIFFUSION APPROXIMATION 


Often in higher dimensional systems, diffusion processes are adopted to analyze complex ge¬ 
netic circuits [7-10]. This section present the derivation to the diffusion approximation of the 
process ( 12 ). 

In the fast-degrading mRNA limit [i.e. 7 —» 00 in (1)], diffusion approximation can be ob¬ 
tained by by performing Kramers-Moyal expansion to (12). The corresponding Fokker-Planck 
approximation reads 


dtPn(x) = -d x 


H \ X )—K - l0X I Pn 


-dt 


H ( x ) 




IO 


IO 


Pn 


Ns is a geometric distribution, and the exact expression of the first two moments are 

E[Ns] = B q 
E [IV*] = Bo (1 + 2B 0 ). 


(17) 


(18a) 

(18b) 


When the bursting number is large B 0 1 (typical biological value 10 1 ~ 10 2 in E. Coli [2, 11]), 
we arrive at the final diffusion approximation of ( 12 ): 


dtPn(x) 


~d x [(bH(x) - 7 0 x)p n ] + ia 2 



TofA 
K ) 



(19) 


We finally remark that in the large population limit, b scales O (K x ). A sensible population 
scaling suggests that near the mean-field fixed points, O ( bH(x )) = O ("lox) = O (if 0 ), which in 
turn indicates that O (H(x)) = O (if). It is clear that the diffusion term can then be simplified 
if we neglect demographic noise due to the protein degradation, when Bq^> 

d t Pn{x) = -d x [(bH(x) - iqx) Pn] + 3 2 [( b 2 H (x)) p„] , (20) 

or equivalently the Ito stochastic differential equation 

dX t = [bH(X t ) - 70 X t ] dt + v 7 2b 2 H(x)dW t , (21) 


where dW t is the Wiener process. 

In the main text, we refer the diffusion process (20) to be the diffusion approximation of 
process (1) in the fast-degrading mRNA limit (7 —» 00 ), with a reflective boundary at x = 0. 

A more phenomenological way is to assert the drift and the diffusion terms of the stochastic 
differential equation to be the mean ( 6 ) and variance ( 6 2 ) of the exponentially distributed burst 
size. It can be shown [12] that this approach corresponds to a constant-burst model. In this 
case, the stochastic differential equation is clearly 

dX t = [bH(X t ) - 70 X t ] dt + y/b 2 H ( x)dW t (22) 


Finally, we remark that (20) is derived from expanding the master equation (12) where effect 
that the bursting noise only enhance the population of the proteins, and ( 21 ) clearly over-estimate 
the noise. On the other hand, the phenomenological approach (22) clearly under-estimate the 
low-density bursts of the exponentially distributed kernel. 
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VI. MEAN SWITCHING TIME OF THE PIECEWISE DETERMINISTIC MARKOV 
PROCESS 


On the domain D := {a; : 0 < x\ < x < x < oo}, the backward equation reads 


1 ^ _ ( -7 + [67 - /<»] d x 
H{x) 


7 \ f x,t) 

~H(x) - f{x)d x ) \T 0 (x,t) 


and in the limit with fast-degrading mRNA 7 — > 00, it converges to 


-1 + bd x 


1 


Ti(x,t) 


H(x) - H(x) - f(x)d x J \ T 0 (x,t) J ’ 
where we denote 70 X by f(x). It is elementary to eliminate the variable Tf and obtain 


d 2 T 0 

dx 2 


+ 


H 1 H d 
f b+xdx 


(f) 


dJ 0 
dx 


Define an auxiliary function M (x) and V (x) 

H{x') 1 H(x') d 


M(x) := J 


1 1 dH' 

~bf ~fH~dx J ■ 


dx', 


(23) 


(24) 


(25) 


/(x') b x' dx' \H{x') J 
y/ x \ _ I ( _I_|_ \ e M ( x ') dx ' 

J ybf(x') f(x')H(x') dx ) 

With the expression H(x) := r 0 + r \ X n /{x n + k n ) and /(x) := 70X, M(x) has a closed form 

1 


(26) 

(27) 


M(x) := log 


Now (25) can be expressed as 


e (x n + k n ) n 'ro 


r 0 


r ix n 
x n -\-k n 


dx 


and the formal solution is 


“ / m / \ dTo \ 

37 e (^777 =“ 


dx 


f 1 


1 dH 


\& 7 ox 70 xH{x) dx 


a M{x) 


T 0 {x)=C 0 + C 1 f e~ M ^' ) dx'+ ( e- M ^ x 'W{x')dx' 

J X\ J X\ 


(28) 


(29) 


(30) 


With two constants of integration Cq and C\. 

Since T 0 (x 1 ) = 0, clearly Co = 0. The second constant C\ can be determined by the second 
boundary condition 


dT 

Ti(x 2 ) = 0 -£=>■ -1 = -H(x 2 )T 1 (x 2 ) - f(x 2 )-^-(x 2 ). 

dx 


After some algebra, we arrive at 
C\ = 


-V{x 2 )e~ M ^) - j ^ 


H{x 2 ) V(x')e M ( x ^dx' — 1 


e -M(x 2 ) + H(x ,) p e -M(x>) dx '' 
f(x 2 ) Jx 1 


(31) 


(32) 


6 























For the mean switching times between the high- and the low-concentration mode, we have to 
impose either X\ —> 0 (for initial state in the low mode) or X 2 —>• oo (for the initial state in the 
high mode). 


A. x\ — v 0 limit 


Note that exp [— M(x)\ has a singularity at x = 0. However, V(x) is a well-behave function 
near x = 0, and we claim 


f e M ^y'>V(y)dx' 
Jo 


< oo. 


To show this, first let 0 < e 1, and 

pX pE pX 

/ e~ M ^V(y)dy = / e~ M ^V(y)dy+ / e~ AI ^V(y)dy. 

Jo Jo J £ 

The second term is bounded. As for the first term, since y < e -C 1, we have 


o~M(y) _ 


e* 


< 


yi 0+1 ( y n + k n ) " 7 o 
e^(r 0 +ri) 1 


ro 


ny 


B i 


n. rx + i ■ ia + i 

fc 7o y to y 7o 


with a constant B\ < oo. Similarly, 


m(v) ^ ^+i ( e " + k n ) ny ° D ia+i 

e < y 7 o ----=: B^y^ 

ro 


with a constant B 2 < 00 . Similarly, V(y) can be bounded: 


V(y) = 


1 


1 dH(z) \ e M(z) dz < o r z ^ dz = _^3 


B-I 10 


b^oz 70 zH{z) dz 
with some B 3 < 00 . Finally, we have 

[ e~ M{v) V(y)dy < = 7 ^ 3-7 [ dy < 00, 
Jo zz + i- Jo 


10 + r 

70 


7o 


which establishes our claim (33). 

Next, we proceed to show the following statements: 


lim 

x\— 




= 0, and lim 


p ~M( x ) dx 

• j X 1 


xi-s-o f x3 e~ M ^dx 

Jx 1 

for any X 3 < X 2 - To show this, again we separate the integral 


= 1 


f e M ( y ^dy= f e M ^dy+ f e M ^dy. 

J X\ J X\ J E 


(33) 


(34) 


(35) 


(36) 


(37) 


(38a) 


(39) 


(40) 
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The second term is again bounded, and for simplicity define 


B 5 := e AI ( y ' ) dy < oo. 


As for the first term, we begin with the lower bound of the integrand: 


As a consequence, 


o-M(x) 


> 


r o 


1 


y~r o + (e" + k n ) n ~<o 


n ■ -^6 r n , ^ • 


Jx 1 r 0 \ _70 £-Y0 


and finally 


lim 

xi—tO 


-M(y) , 


< lim 


ro 


r 0 

X^o 


^670 x 7 o + e^o 


= 0 . 


( 41 ) 


(42) 


(43) 


(44) 


Similarly, it is straightforward to apply the L’Hopital’s law to show 

f x2 e~ M ^dx 

lim -——— = 1. 

antO f e~ M ^dx 

JX 1 


(45) 


To sum up, upon taking the limit x\ —> 0, the solution (30)—the mean switching time if the 
system starts with the low-protein-abundance mode—can be expressed as 


with 


T 0 (x)= lim To(x) = [ e~ M{x ' V(z')dz' + T o (0) 
zi-s-0 J 0 


T 0 (0) = lim Ci / e" M(:r 'W 
Jxi 


= lim 

xi—>-0 


-V{x,)e- M ^ - 


H{x 2 ) y (x')e~ M ^ x '^dx' - 1 


e -M(* 2 ) + jU£2> r*2 e -M(x') dx , 
_ fj; j 1 ; 1 _ 

f x e- M ^ , '>dx' 

Jx 1 


~V(x 2 )e ~ M ^ 


fffca) / X2 V(x')e- M ^dx' - 1 


H(x 2 ) 
f( x 2 ) 


1 


ff(cc 2 ) L 


1 - x 2 V{x 2 )e~ M ^ ) 


rx 2 

/ e~ M ^ x '' ) V(x')dx' 

Jo 


(46) 

(47a) 

(47b) 

(47c) 

(47d) 























B. X 2 —>• oo limit 


Note that exp [—M{x)\ has a singularity at x —>• oo. On the other hand, V{x) is again well¬ 
behaving. Rewrite 


C X = 


-V(x 2 ) - j^j eM( - X2) H(x 2 ) f X2 V{x')e M ^dx' - 1 


l + e M( X2 )m?2} r*2 M ( X ') dx , 

f(X 2 ) Jx 1 

and we aim to show that hin^^oo C\ = — liim^oo V(x) < oo. 
First, we show that for a given and strictly positive x\, 

rX 2 

lim e M(x2) / e~ M{y) dy = 0. 

J X\ 


Clearly from the exact expression of M(x 2 ), we have 


^-+1 


0 M(x 2 ) 


px 2 
J X\ 


xJ 0 {x% + k n )^ 


e ~ M (v)dy = I e - {x2 ~ v) 


r o + - 


y 7 ° +1 ( y n + k n ) "^0 \ iy 

r 0 ' v n +, 


-dy 


< 1 i + iiWi + 

ro 


k n \ ^ 

*^2 / Jx\ 


y' n +k Tl 

IQ. + IL +1 
3 ; \ 70 70 


(48) 


(49) 


(50) 


dy. (51) 


For x 2 ~3> Xi, it is elementary to show that the integrand has a single maxima at y* := ro /70 + 
fi /70 + 1 < 00 . As a consequence, 


0 M(x 2 ) 


rx 2 

/ e~ M ^dy 

J X\ 


< 1 


Tj_ 

ro 


r 1 


1 + 


\ n 70 


,y*-x 2 


^ + — + 1 /.X9 

^2 \ 70 70 / 2 


= I 1 + — I I 1 + — 

7-0 


\ n 70 


V* 


X 2 \ TO 70 


dy 


(52a) 


2/* 


{x 2 -xi), (52b) 


and apparently, 


lim e M(x2) / e~ M{v) dy = 0 . 

x 2 —^OO ' 


(53) 


Note that V(x) is a strictly decreasing and well-behaved function, so 0 > V(x) > Vmin := 
lim :c _ ) . 00 V(x) for x\ < x < 00 . As a consequence, 


lim 

x 2 —too 


0 M(x 2 ) 


px 2 /*x 2 

/ V{y)e~ M{ - v) dy < |R min | lim e^ 2 } / e~ M{v) dy = 0. 
Jxi X2 ^°° Jx 1 


The above Eqs. (53) and (54) suggest that 


lim Cl -7 |Vmin| , 
x 2 —too 


(54) 


(55) 


and the final solution of the mean switching time if the system begin with a high-protein- 
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abundance state is 


To Or) 


IK 


min | / & 

J Xi 


-M(x 


^ dx' + / 

J X-i 


e - M{x ') V {x’) d x'. 


(56) 


VII. MEAN SWITCHING TIME OF THE DIFFUSION APPROXIMATION 


The mean first-exit time of general one-dimensional diffusion process was extensively studied 
in the classic book by van Kampen [4]. For the reference of the reader, in this section we briefly 
present the key equations used to evaluate the mean switching time of the process ( 21 ) 

Given any transition boundary x c , we define the domains of interest 0^ := {x : x < xc } and 
fin '■= {x : x > x c }, respectively the low and high protein abundance states. For the domain 
Ql, we impose a reflective boundary at x = 0. We are interested in the mean of the first times 
when the process X t leaving the domain H £ {H l, ^h}- 

Denote the random first exit time by t (x) := inf {t, : X t <£. D|Xo = x}. The mean first exit 
time T(x) := E [r (x)] satisfies the following adjoint equation 

- 1 = + D(x)^T(x). (57) 

Here, we define the drift v(x) and the diffusion D(x) according to (21): 

v(x) := bH(x) — 70 X, 

D(x) := b 2 H (x). 

The general solution of (57) is 

T (x) = C 0 + C\G(x) - H(x) (60) 

where Co anc C\ are constant from the integrations, and the auxiliary functions are defined to 
be 


(58) 

(59) 


x o) 


G(x;x 0 ) 


K(x\xo) 


H(x; xq) 



f e~^ v) K ( y) dy. 

J Xn 


(61) 

(62) 

(63) 

(64) 


The boundary conditions of (60) depends on which domain the initial state is on. On the one 
hand, when the state start from a low protein-abundance state, the boundary conditions are 

jm 

—— (0) = 0, and T(x c ) = 0 (65) 

dx 

which result in the mean switching time from low to high protein-abundance to high protein- 
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abundance state 


T(x < x c ) 


lim H (xi',x c ) 

X\ —^OO 


G(x;x 0 ) - H(x-x o). 


( 66 ) 


On the other hand, when the state begin with a high protein-abundance state, the boundary 
conditions are 

dT 

lim -—(a;) = 0, and T(x c ) = 0, (67) 

x —>oo dx 

and the mean switching time from high to low protein-abundance state is 

T{x<x c )=H(x c ;0)-H(x]0). (68) 
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